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Abstract 

The usefulness of piecewise polynomials with C 1 and 
C 2 derivative continuity for response surface 
construction method is examined. A Moving Least 
Squares (MLS) method is developed and compared 
with four other interpolation methods, including 
kriging. First the selected methods arc applied and 
compared with one another in a two-design variables 
problem with a known theoretical response function. 
Next the methods are tested in a four-design variables 
problem from a re liability- based design application. In 
general the piecewise polynomial with higher order 
derivative continuity methods produce less error in the 
response prediction. The MLS method was found to be 
superior for response surface construction among the 
methods evaluated. 

Introduction 

Structural reliability engineering analysis involves 
determination of probability of structural failure taking 
into account the uncertainty in the geometric 

parameters, material properties and loading conditions'. 
The uncertain quantities are treated as random variables 
with known probability distributions in probabilistic 
analysis. In these analyses, for each set of random 
variables a Monte Carlo simulation is performed to 

determine the probability of failure of the structure . 
Monte Carlo simulations require large number of 
simulations and hence, are computationally expensive 
and also require large amounts of analyses time. In 
order to alleviate some of the problems associated with 


Monte Carlo simulation, approximate methods such as 
First Order Reliability Method (FORM) and Second 
Order Reliability Method (SORM) were developed to 

minimize the number of simulations . Even with these 
approximate methods, the computational effort required 
to perform a structural reliability analysis can be very 
high. Thus, there is a need to develop methods that are 
accurate and yet minimize the computational time. 

Response surface functions are often used as simple 
and inexpensive replacements for computationally 
expensive structural analyses in reliability methods. In 
the response surface method, a surface is fit to 
interpolate data that are generated by performing 
structural simulations at selected points in the 
parameter space. The response surface is then used to 
evaluate the structural response at other points in 

parameter space. For example, the polynomial 

3 

regression response surface method is widely used in 
structural optimization. The number of data points 
necessary to fully sample the parameter space increases 
exponentially as the number of random variables or 
design variables increases. In theory it is possible to 
find an interpolating polynomial of sufficiently high 
order to pass through all the data points in the 
parameter space. However, for such interpolating 
polynomials as the number of data points increases, the 
surface generated tends to exhibit oscillations or 
wiggliness” between the data points. For this reason the 
polynomials are generally limited to first and second 
order to represent the response surface, and Least 
Squares regression is used to best fit the data points 
Many other methods use piecewise polynomials instead 
of a single global polynomial to represent the entire 
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parameter space. The piecewise polynomials arc 
generally limited to C° continuous surfaces. In C° 
continuous surfaces, the function values are continuous 
everywhere (across piecewise polynomials), hut first 
derivative (slope) continuity and higher order 
derivative continuity are not guaranteed. Recently 
Romero el al. 4,5 implemented a Progressive-Lattice- 
Sampling fPLSl method based on finite element 
interpolation ( C° globally continuous) to construct the 
response surfaces. However, response surface 
construction methods using piecewise polynomials with 
C 1 continuity (function values and slopes are 
continuous) and with C 2 continuity (function values, 
slopes and second derivatives are continuous) are rare 
in the literature. 

The purpose of the present paper is to study the 
usefulness of piecewise polynomials (with higher order 
derivative continuity) on response surface generation 
and its application to reliability engineering. The 
Minimum Norm Network (MNN) developed in 
references 6-10 adopted to construct a response surface 
with C 1 and C 2 continuity in two design variables. A 
new enhanced C 2 MNN method is developed to 
reproduce polynomials up to third degrees exactly. 
Also a Moving Least Squares (MLS) method is 
developed to reproduce C and C' continuous 
response surfaces for arbitrary number of design 
variables. The selected interpolation methods are tested 
in a two variables problem and results are compared. 
Finally a four-variable example from a reliability 
application 11 is presented to demonstrate the 
effectiveness of the MLS method. The MLS method is 
compared with least squares polynomial and the kriging 
methods on this problem. 


Interpolation Methods 

The accuracy of the response surface in representing 
the behavior of the actual system largely depends upon 
the interpolating method used for its generation. Brief 
introductions to the five interpolating methods used in 
the current analysis are presented in this section. 


example, a quadratic polynomial with NDV design 
variables has the form 

a v>r A7)r 

y(*) = * o + !*,■*,■+ I MV V / «> 

7 = 1 7 = 1,/=/ 

Where j> is the approximated value of the target 

function at the point in the parameter space having 
coordinates (X } , X 2 X NDV ) , and a {] , a ■ , and h tJ are 

the unknown constant coefficients. The unknown 
coefficients are determined by a regression procedure. 
Most commonly, the method of least squares is used to 
determine the coefficients that minimize the error of the 
approximation at the sampling points . Since a single 
polynomial is used to represent the entire parametric 
space, the method is termed here as the Global Least 
Squares (GLS) method. In the present study, the GLS 
method is limited to the quadratic polynomial given by 
Equation (1). 

2. Kri gi n g 

Kriging is an interpolation method that originated in the 
geostatistics community. Kriging uses the properties of 
the spatial correlation among the data samples. In 
arriving at an interpolated value at some point in the 
parameter space, kriging more heavily weights data 
samples that are “nearby” rather than giving all data 
samples equal weight. This is achieved by setting mean 
residual error to zero and also minimizing the variance 
of the errors. The final equations for kriging are given 
below from reference 12 for N sampling points 
and NDV design variables: 

The estimated value of y in kriging is obtained from 

y = P + r‘ R-'iy-fif) (2) 

where Y is the column vector of known function 
values at the N sampling points, /? is a constant to be 
determined, R is correlation matrix obtained for an i th 
row and /^column from the correlation function as 


R(X\X J ) = exp 


NDV 

-e i x k 

k~\ 



(3) 


L Global Least Squares Method (GLS) 

The GLS methods are generally known as polynomial 
regression methods and are widely used in the literature 
1 3). The GLS methods are used to create response 
surface functions from a set of sampling points. For 


r is the column vector of length N obtained 
from 


r‘ = |/?(A',A’ 1 ), R(X,X 2 ) R{X,X N )} (4) 
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and / is vector of length N , with all elements in the 
vector set to unity as 

/' if (5) 

The unknown /? in Equation (2) can he obtained from 

fi=(f r R~ l fy l f r R-'y ( 6 ) 

The Maximum Likelihood Estimate (MLE) for the 
unknown quantity 0 in Equations (3) is obtained from 
a one-dimensional maximization problem defined by 


MLE= max 



0 < $ < oo (7) 


derivatives at the vertices. The unknown 
partial derivatives are determined by enforcing 
continuity of all edges meeting at a vertex by 
minimizing a functional (second order for 2 
parameters). The union of all the edges of all 
the triangles forms the MNN network. At the 
end of the second step, function values and its 
derivatives are known at the vertices of the 
triangles in the network. 

3. The interpolation is extended inside the 
triangles in terms of the function values and its 
derivatives at the vertices of the triangle. The 
interpolation inside the triangles is achieved 
using triangular blending functions forming a 
piecewise continuous surface. The blending 
function is constructed in such way that it 
produces C ! and C 2 continuous surface with 
its neighboring triangles. 


Where 

^2 = (r-i/f/r'Cr-i/ ) ^ 


Estimation of 6 in the one dimensional optimization 
problem is the critical step in the kriging method. The 
kriging method used in this study produces a C 2 
continuous interpolating function over the entire 
parameter space. 

3. Minimum Norm Network (MNN): 

The MNN method is a piecewise interpolation scheme 
that can produce C 1 and C 2 continuous surfaces. 
Complete details about the method can be found in 
references 6-10, only a brief description of the method 
is presented here. The MNN method requires the 
function values at the arbitrary sampling points. The 
procedure to construct a MNN can be described in the 
following three steps: 

1. The given sampling points are triangulated 
such that each point lies at a vertex of the 
triangles. 

2. The derivatives are generated at the vertices of 
the triangles from known function values. On 
each edge of the triangles an interpolating 
polynomial is fit in terms of the known 
function values and unknown partial 


A response surface with C ] and C 2 continuity is 
produced using the method developed in references 6-9 
and used for the interface element development in 
reference 10. The MNN from reference 8 produces a 
C 2 continuous surface in the limit. Hence, a new 
“enhanced C 2 ” method is developed here to produce 
C continuity (partial derivates up to third degrees 
continuous) in step 2 of MNN network and uses the 
same blending function as that of C 2 in step 3. 

At present the MNN method is available only for two 
design variables (dimensions). Extending the MNN to 
more than two variables is difficult. However, these 
methods are included in this study, since, they are the 
only methods that provide piecewise polynomial 
interpolation. Also these methods can produce a 
C 2 continuous surface that passes through all the 
sampling points for arbitrarily oriented data points (i.e., 
not in a structured grid). 


4. Moving Least Squares Method (MLS): 

The Moving Least Squares (MLS) method is widely 

13,14 

used in meshless methods . Recently the MLS 
method has been successfully applied for response 
surface generation in the context of optimization in 
reference 15. A MLS method is developed here for an 
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arbitrary number of design variables. The method is 
briefly discussed here: 

The MLS approximation for the estimated value y(x) 
can be written as 

y{x\=p r {x)a m {X) (9) 

where p 1 {x}=[p l {X},p 2 {X},...p m {X}] is a 
polynomial basis function of order m used in the MLS 
interpolation. a m (X) is a vector containing coefficients 
a f (X\ j = 1,2,. ..,iw , which are functions of the spatial 
coordinates. For example, for two design variables, 


p ! (x} = [UX u X 2 \ Linear basis function; m- 3 
p‘ \x}= [lX l ,X 2 ,X^,X i X 2 ,Xl\ ( 10 ) 

Quadratic basis function; m - 6 


The unknown coefficients a m can be determined using 
the weighted least squares error norm J(X) at the N 
sampling points 

J(X)^ Z^(X)[p r (X,)a m (X)-Y\ () n 

= [Pa n AX)-Y] w[Pa m {X)-Y] 

where w,(A) is weight function associated with node 
i , whose value is nonzero only in the support or 
influence domain of the node X, ( usually a sphere of 
radius R , ). The matrices P and W are defined as 


P = < 


p'\X i) 

p'\x 2 ) 


p‘\x N ) 


.Vxm 


( 12 ) 


H’i(A') 


W = 


0 


0 


W N (X) 


JNxN' diagonal matrix 


and 


(13) 


04) 


Minimizing the norm J(X) in Equation (10) with 
respect to a m (X) leads to the following linear relation 
between a„,(X) and y 

A(X) a m (X)= B( X ) Y (15) 

where the matrices A(X ) and 5(A) are defined by 


A(X)= P r WP= iw,{X) p(X,) p' iX,) (16) 

i~ 1 

B(X) = P'W 

= [w,(^) / 7(J)r 1 ) w 2 {X)p{X 2 ) - u, v (A')p(A-, v )] 

(17) 

The unknown coefficients a m (X) can be obtained by 
solving Equation (15), which results in 

a n ,(X) = A\X)B{X)Y (18) 

Substituting the unknown coefficient from Equation 

(18) into the Equation (9) leads the MLS interpolation 
of the y as 

y~p'XX)A-\X)B{X)Y (19) 

The MLS approximation given in Equation (18) is well 
defined only when the matrix A is non-singular. This is 
true only if there are at least n sampling points in the 
influence domain of a node X , such that n>m. For 

example, for a one-dimensional case with a linear basis 
function {m = 2) , the value of n should be > 2 . For a 

quadratic basis function in a two-dimensional case the 
value of n should be > 6. 

Except for the weight function vr,( X) , all other 

quantities in the MLS approximation are well defined. 

As already mentioned, the weight function is non-zero 
only in the influence domain of a node i , and equal to 
zero outside the influence domain. In the present study, 
the influence domain is assumed to be of a sphere with 
radius / , . The radius l, must be large enough to 

contain at least m nodes in each direction of the 
parametric space. The weight function is selected such 
that its value goes from unity at the center of the 
influence domain to zero at the boundary and outside 
the influence domain. This property of the weight 
function makes MLS a local approximation compared 
to the GLS approximation traditionally used to 
represent the entire domain by a single function. It may 
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be noted that in the MLS method for every new 
interpolation point (y) Equation (18) is formed and 
solved. 

In this paper, three spline functions with cVC 2 , and 
C 3 continuity are used as weight function 

For C' : 


w,(x) = 


\-3p?+2tf 


0<p,<\ 

P,>\ 


( 20 ) 


For C~ 


w,(x) = 


I - 1 Op, + 1 5/7, 4 - 6/7, 5 


0<A<l 
n, >*i 


( 21 ) 


For C 


•i3 


”’,(*) = 


1 -35 p, 4 +84p, 5 -70/0," ■* 2 0/9, 7 


0 < Pj > 1 

Pt > >, 

(22) 


\J\ j — J\ \ 

where p, = J L is the normalized distance, from 

i, 

the center of the influence domain ( X t ) and a general 
point X . 

The smoothness of MLS approximation is controlled by 
both the weight and basis functions. The precision 
(continuity) of MLS interpolation will be equal to the 
minimum precision of the weight and basis function. 

5. Piecewise Finite-Element (FE) Interpolation 

Piecewise finite element (FE) interpolation was 

. . 4 

originally used in conjunction with Progressive Lattice 
Sampling (PLS). PLS is discussed later in this paper. 
The arrangement of samples in each PLS Level allows 
the parameter space to be subdivided into a regular 
pattern of adjacent polygons, which for, two- 
dimensions, results in triangular and quadrilateral finite 
elements. For the 2-D parameter space, a combination 
of 3-, 4-, and 6-node triangles and 9-node quadrilaterals 
exist at the various PLS Levels. Low-order polynomial 


interpolation is applied over each clement. The 
collection of all the elements together creates a locally 
compliant C° continuous function over the parameter 
space. A detailed discussion of the implementation of 
the FE interpolation technique used here is presented in 
reference 4. 


X, 


X, 


Level- 1 (3 samples) 



Level-3 and 4 (+4=9) 



Level-6 (+12-25) 


4 


' • ' 

Level-2 (+2=5 samples) 


• 

• 

• 

• 

Level-5 (+4=13 

f+J + 

4 + 1 

tffL 

it" t / 

+ t+l, 

W t* 

+.+ 1 

i • • 

• • i 


Level-7 (+16=41) 


Figure 1. Progressive Lattice Sampling Points 


Progressive Lattice Sampling (PLS) Experimental 
Desig n: 

The selection of sampling points plays a major role in 
the accuracy of a response surface. There are many 
schemes available in the literature (see reference 3). 
Romero et al. 4 used Progressive Lattice Sampling 
(PLS) incremental experimental design sequence as 
shown in Figure 1 for two variables X } and X 2 . In 

this example, the square represents the parameter space 
of X , and X 2 . Level 1 of the design consists of three 

samples, with one sample in the center of the parameter 
space and two other samples along the boundary. For 
n parameters. Level I requires n + 1 samples. Level 2 
adds n samples to complete a 2^+1 “simple- 
quadratic" layout. Level 3 adds a 2 ” factorial design. 
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Level 4 adds a Box-Behnken design to complete an 
over all 3" full factorial design. (In 2-D, Levels 3 and 
4 have the same layout.) Level 5 adds a sub-scaled 2 " 
factorial design as shown in the figure. Level 6 adds the 
appropriate samples to complete a 5" full factorial 
design. Level 7 adds a sub-scaled 4" full-factorial 
design in the interior of the parameter space as shown. 
The strength of PLS is that it provides an efficient way 
to add sample sites that leverage previous samples so 
that uniform distribution of the samples over the 
parameter space is maintained. The PLS design points 
will be used as sampling points in all the fitting and 
interpolation schemes in the present study. 

A pplication Problem 
Two-Variables Problem: 

First the described interpolation methods are applied to 
a two-design variables problem selected from reference 
4. Target response function for the two-variables is 
shown in Figure 2 



Figure 2. Target Response Function for Two- 
Variables Problem 


To examine the fitting performance (within the PLS 
frame work) of the various response surface 
construction methods, a global measure of average 
error is defined as follows: 

■v 

X | {exact ), - ( predicted ), | 

A verage Error = — ( 24 ) 

N 

Where “exact” in the summation comes from the 
evaluation of the exact function. The predicted values 
in the summation come from the response surface 
approximation at V interpolated points For this 
example Vis set to equal to 441 and selected from 
equally spaced points on a 21x21 square grid overlaid 
on the domain. Earlier experience in reference 5 
indicates the 21x21 grid appears to be sufficiently 
dense to achieve adequate representation of the target 
and approximate functions. 

Four levels with 9, 13. 25 and 41 points were selected 
for comparison. For this example problem, the average 
errors shown in the following figures for kriging and 
finite element interpolation are taken from reference 5. 



Figure 3. Average Errors in Two- Variables Problem 
in Level-4 with 9 Sampling Points 


This response function is defined as: 
response( , X ? ) 


0.8r + 0.35S7/7 


2.' 
in — 


Am 


[l ,5.S7w(l .3(9)] 


l V2 

on the domain 0 < < 1, 0 < X 2 < 1 


JX? + X} , 6 = arc tan 

f* 2 ! 

V 1 - 

U.J 


Figures 3 to 6 show the average errors from levels 4 to 
7. Since 9 points are not enough to fit a quadratic basis 
function in MLS interpolation, a linear basis function is 
used. Similarly the 9 points are not sufficient to fit the 
Enhanced C : -MNN and hence it is not shown in Figure 
(23) 3. Finite element results for level-4 were not available 

and are not shown in Figure 3. 


Exact data values of this function are obtained at the 
PLS sampling points shown in Figure I. The response 
surface is generated for each of the various PLS levels 
in Figure 1. The response surface is then used to 
interpolate the value at any other point. 
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Figure 4. Average Errors in Two-Variables Problem 
in Level-5 with 13 Sampling Points 



Figure 5. Average Errors in Two- Variables Problem 
in Level-6 with 25 Sampling Points 


0.25 



Figure 6. Average Errors in Two-Variables Problem 
in Level-7 with 41 Sampling Points 


From Figures 3 to 6, the following observations can be 
made: 


2. At lower number of sample points, level-4, all 
the methods produce almost the same average 
error, 

3. In all the levels the GLS method consistently 
performs poorly. This implies that the GLS 
method should be avoided for response surface 
generation, 

4. All the piecewise interpolation methods 
including the finite element interpolation 
methods perform better at all levels, 

5. The locally weighted methods, kriging and 
MLS, perform as well as the piecewise 
interpolation method like the MNN method, 

6. The C'-MNN method performed as well as C - 
MNN method at all levels, 

7. All the methods with higher order derivative 
continuity produce less error at all levels, and 

8. The kriging and Enhanced C - MNN methods 
produce the least error at level 7, closely 
followed by MLS method. 

Since it is easy to extend the kriging and MLS methods 
for arbitrary number of variables, the two methods are 
compared separately in Figure 7. It can be concluded 
from Figure 7, that MLS consistently produced 
marginally smaller errors than kriging at all levels 
except level-7 with 41 sampling points. However, since 
the difference in error between the two methods is 
small, there is no clear advantage of one method over 
the other. 

0.25 

u 0.2 

£ 0.15 

CL) 

8 o.i 

> 

< 0.05 
0 

0 


Figure 7. Kriging and MLS Methods Comparison 
for Two-Variables Problem 



Number of Sampling Points 


1. All the interpolation methods show reduced 
average errors as the number of sampling 
points increases, 


Four- Variables Problem: 

The next example problem is taken from reliability- 
based design of a metal, plate- 1 ike wing to meet 
strength and flutter requirements". The selected plate- 
like wing configuration is shown in figure 9. 
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The dimensions used for wingspan (L ) , wing root 
chord (C }i ) , tip chord (C, )* and sweep of the leading 
edge (O) are also show in Figure 9. The modulus of 
elasticity is 1 0 x 1 0 6 psi and Poisson's ratio is 0.30. 

The wing is clamped at the root and subjected to a 
uniform pressure of 1 psi . 


I 4 



where -1<£<1, and -!<;;<] 


Element 1 



Element 162 


Figure 9. Finite Clement Model for Stress Prediction 


▼ x 

Figure 8. Dimension of Metal Plate-Like Wing 

The thickness distribution along the span of the wing is 
assumed as bi-linear and can be defined in terms of the 
thicknesses of the corner nodes I to 4 (see Figure 8) as 


= c, + c 2 4 + ctf + c^T] 

( 25 ) 

where 


(6 +h +6 + 6 ) 

C|_ 4 

(26) 

(-/, +t 2 +/ 3 “/ 4 ) 
4 

(27) 

(” 6 “* 2 +* 3 + * 4 ) 

C-L ~ 

4 

(28) 


c 4 


(/} - t 2 + / j - *4 ) 
4 


(29) 


The equation relating the (x, y) and wing 

coordinates (see Figure 8) can be written as 


2L(x~ y tan 0) ^ 

C r L-(C r -C f )y~ 



(30) 

(31) 


The four corner node thicknesses (/, to t A ) are the 

design variables. Each thickness is allowed to vary 
between 0. 1 5/>? and0.4/>? . The sampling points are 

generated using the PLS scheme for level-5 and level-7. 
There are 97 PLS sampling points for level-5 and 88 1 
PLS sampling points for level-7. In order to predict the 
stress distribution as function of four design variables 
(/] to t 4 ), the plate is divided into 162 quadrilateral 

finite elements as shown in Figure 9. Finite element 
analyses with 162 quadrilateral elements are used to 
obtain the stresses at the centroids of each element. For 
example, for the 881 sampling points in level-7, 881 
finite element analyses are performed. These 881 
cetroidal stresses for an element are used to construct 
the responses for that element. Hence, total of 162 
response surfaces are constructed, one for each 
element. The NASTRAN structural analysis code with 
8-node quadrilateral elements is used for the finite 
element analyses. 

The 162 response surfaces (one for each finite element) 
are generated using GLS, kriging, and MLS methods in 
terms of the design variables t } to t x . The fitting of the 

response surface with GLS and MLS is straightforward, 
but, the kriging method requires some clarification. As 
previously mentioned, the estimation of 6 in Equation 
(4) using one-dimensional optimization is the critical 
step in kriging. The variation of Maximum Likelihood 
Estimate (MLE) with 0 is shown in Figure 10. 
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-9600 

Figure 10. Variation pf Maximum Likelihood 

Estimate (MLE) with Free Parameter 0 

It can be seen from Figure 10 that there is no absolute 
maximum value of MLE. The MLE reaches the 
maximum value at the default minimum value of 0 = 0 
degrees. The matrix R in Equation (3) becomes 
singular for 0 = 0 degrees. For the current problem 
with 881 sample points, the determinant of the matrix 
R becomes numerically nearly /ero for 0< 3. The 
selection of 6 here then becomes highly subjective. 
The value of 0 is set to 0 = 1 and 0 = 5 for level-5 
and level-7 respectively. These 0 values are selected, 
since they produced accurate results. 
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Error 



Figure 11. Average Error Distribution for GLS 
Method Over the Wing Span 


Error 



To examine the fitting performance, 2500 randomly 
selected thickness sets are used in Equation (23) to 
obtain the average error. The average errors are 
calculated for alt the 162 elements. The stress values 
obtained from the 2500 NASTRAN analyses are used 
as ‘exact' in Equation (23). The average error for each 
element is normalized by the average stress for that 
element using 


25’K) 

X | exact stress\ 

Mean Stress for an element = — (32) 

2500 

The percent error is calculated as 


%Error for a n element . * l0 0 Wl 

Mean Stress 


Figure 12. Average Error Distribution for Kriging 


Error 



0 245 I 

0.23°I 
0.2 1 5§f 
0.200 


0 1 85 
0 . 1 70 

0 155: 
0.140-’ 
0 125: 

o no 1 

0.0951 
0.0801 
0 0651 
0 0501 
0.034* 


Method Over the Wing Span 


Figure 13. Average Error Distribution for MLS 
Method Over the Wing Span 
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The distributions of average errors over the wingspan 
are shown for GLS, kriging and MLS methods in 
Figures 11 to 13. The relative differences in average 
error for the three methods are almost the same in each 
element. As a representative example, the percent error 
in Element 1 located near the root of the wing (see 
Figure 9) is compared for the three interpolation 
methods in Figures 14 and 15 for levels 5 and 7. For 
level-5 with 91 points, the GLS produced 8.65 percent 
error compared to 1.8 percent for kriging and 1.5 
percent error for MLS. The corresponding values for 
level-7 with 881 points are GLS=5.8 percent error, 
kriging=1.7 percent error and MLS=0.58 percent error. 
The kriging and MLS methods produce errors that are 
an order of magnitude less than the GLS method. 
Numerically the MLS produced the least error at both 
levels. However, between kriging and MLS, it is 
difficult to choose one over the other. Further study is 
needed to draw definite conclusions. 

Discussion 

In the above two examples it was shown that the 
performance of the kriging and MLS methods are 
nearly the same. However, in kriging there is a need to 
estimate the free parameter 0 through optimization. In 
the MLS method there is no parameter to evaluate, 
except to define the influence radius /, . The free 

parameter l } is easy to select from the requirement of 

the number of points to make the matrix A in Equation 
(15) non-singular. It should be noted that the kriging 
and the MLS methods are called locally weighted 
interpolation rather than piecewise interpolation. It was 
observed in the two variables example that the 
piecewise polynomials produced small error at all the 
levels. But these methods are available only for two 
variables. Hence, piecewise polynomial methods for 
arbitrary number of dimensions would need to be 
developed. 

Concluding Remarks 

The usefulness of Piecewise Polynomials with higher 
order derivative continuity methods are studied and 
evaluated for response surface construction. The 
piecewise polynomials with higher order derivative 
continuity (C'and C~ continuity) methods produce 
fewer errors for a given set of sample points than the 
global least squares methods. The kriging and MLS 
methods performed equally well and are easy to extend 
to arbitrary numbers of design variables. However, it is 


easier to fix the free parameters in the MLS method 
than in the kriging method. There is a need to develop 
or examine piecewise polynomial methods that can be 
easily extended to arbitrary number of design variables. 
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Figure 14 . Average Errors: Four- Variables 
Problem With 97 Sampling Points 
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Figure 15. Average Errors: Four- Variables 

Problem With 881 Sampling Points 
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